setwd("C:/Documents and Settings/wjm02/Documents/brca_metabric/")
##### data format ##### -----
library(dplyr)

exprset = read.table("data_expression_median.txt",header=T,sep = "\t")

matrix = as.matrix(exprset[,3:ncol(exprset)])
rownames(matrix) = exprset[,1]
tmp <- by(matrix,
          rownames(matrix),
          function(x)x[which.max(rowMeans(x)),])
expr <- do.call(rbind,tmp)
write.table(expr,"expr.txt",sep="\t")

meta_patient = read.table("data_clinical_patient.txt",header=T,sep="\t")
meta_sample = read.table("data_clinical_sample.txt",header=T,sep="\t")
tnbc_sample = meta_sample[meta_sample$ER_STATUS=="Negative"&meta_sample$PR_STATUS=="Negative"&meta_sample$HER2_STATUS=="Negative",]
tnbc_patient = left_join(tnbc_sample,meta_patient,by="PATIENT_ID")
tnbc_patient = tnbc_patient[tnbc_patient$ER_IHC!="Positve",]
tnbc_patient = tnbc_patient[tnbc_patient$VITAL_STATUS!="Died of Other Causes",]
tnbc_patient = tnbc_patient[,c(1,2,7,11,25,26,29)]
tnbc_patient$SAMPLE_ID = gsub("\\-", ".", tnbc_patient$SAMPLE_ID, ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
tnbc_patient$OS_STATUS = gsub("1\\:DECEASED", "DECREASED", tnbc_patient$OS_STATUS, ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
tnbc_patient$OS_STATUS = gsub("0\\:LIVING", "LIVING", tnbc_patient$OS_STATUS, ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
write.table(tnbc_patient,"tnbc_patient.txt",sep="\t")

##### TNBC ##### -----
rm(list=ls())
expr = read.table("expr.txt",sep="\t",header=T)
tnbc_patient = read.table("tnbc_patient.txt",header=T,sep="\t")
tnbc_sample = data.frame(ID=(tnbc_patient$SAMPLE_ID))
expr_sample = data.frame(ID=(colnames(expr)))
tnbc_sample = merge(tnbc_sample,expr_sample,by="ID")
tnbc_sample = tnbc_sample$ID

tnbc_expr = expr[,tnbc_sample]
tnbc_zscore = t(scale(t(tnbc_expr)))

tnbc_sample = data.frame(SAMPLE_ID=tnbc_sample)
tnbc_patient = left_join(tnbc_sample,tnbc_patient,by="SAMPLE_ID")

write.table(tnbc_patient,"tnbc_patient.txt",sep="\t")
write.table(tnbc_expr,"tnbc_expr.txt",sep="\t")
write.table(tnbc_zscore,"tnbc_zscore.txt",sep="\t")

##### survival ##### -----
rm(list=ls())

tnbc_sample = read.table("tnbc_patient.txt",header=T,sep="\t")
tnbc_zscore = read.table("tnbc_zscore.txt",header=T,sep="\t")

tnbc_zscore = as.matrix(tnbc_zscore)

library(GSVA)

IFNA_gene = data.frame(gene=c("HLA-G","PTPN2","IRF9","LGALS3BP","NUP210","EIF2AK2","KPNA7","CXCL10","NLRC5","UBE2L6",
                              "B2M","IFI6","IL-15","IRF1","STAT1","C1S",
                              "GBP6","GMPR","USP18","LAMP3"))
IFNA_score = gsva(tnbc_zscore,IFNA_gene)

library(survival)
library(survminer)
rownames(tnbc_sample) = tnbc_sample$SAMPLE_ID

high_score = as.data.frame(IFNA_score[,IFNA_score["gene",]>0])
colnames(high_score) = "score"
IFNA_high = tnbc_sample[rownames(high_score),]
IFNA_high$group = "high"
IFNA_high = cbind(IFNA_high,high_score)
IFNA_high = IFNA_high[order(IFNA_high$score),]


low_score = as.data.frame(IFNA_score[,IFNA_score["gene",]<0])
colnames(low_score) = "score"
IFNA_low = tnbc_sample[rownames(low_score),]
IFNA_low$group = "low"
IFNA_low = cbind(IFNA_low,low_score)
IFNA_low = IFNA_low[order(IFNA_low$score),]


dat = rbind(IFNA_high,IFNA_low)
write.csv(dat,"dat.csv")

my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECREASED')
kmfit1 <- survfit(my.surv~group,data = dat)
curve = ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE,legend.title="",xlab = "Months") +
  ggtitle("IFNA score survival curve")

pdf("./metabric IFNA score.pdf", onefile = FALSE)
print(curve)
dev.off()
curve




